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ABSTRACT 

We use one-dimensional two-zone time-dependent accretion disk models to study the long-term 
evolution of protostellar disks subject to mass addition from the collapse of a rotating cloud core. 
Our model consists of a constant surface density magnetically coupled active layer, with transport 
and dissipation in inactive regions only via gravitational instability. We start our simulations after 
a central protostar has formed, containing ~ 10% of the mass of the protostellar cloud. Subsequent 
evolution depends on the angular momentum of the accreting envelope. We find that disk accretion 
matches the infall rate early in the disk evolution because much of the inner disk is hot enough to 
couple to the magnetic field. Later infall reaches the disk beyond ~10 AU, and the disk undergoes 
outbursts of accretion in FU Ori-like events as described in Zhu et al. 2009c. If the initial cloud core is 
moderately rotating most of the central star's mass is built up by these outburst events. Our results 
suggest that the protostellar "luminosity problem" is eased by accretion during these FU Ori-like 
outbursts. After infall stops the disk enters the T Tauri phase. An outer, viscously evolving disk has 
structure that is in reasonable agreement with recent submillimeter studies and its surface density 
evolves from S oc Rr 1 to Br 1 ' 5 . An inner, massive belt of material- the "dead zone" - would not have 
been observed yet but should be seen in future high angular resolution observations by EVLA and 
ALMA. This high surface density belt is a generic consequence of low angular momentum transport 
efficiency at radii where the disk is magnetically decoupled, and would strongly affect planet formation 
and migration. 

Subject headings: accretion disks, stars: formation, stars: pre-main sequence 



1. INTRODUCTION 

In the picture of low-mass star formation from large, 
cold protostellar clouds, any small initial rotation of the 
cloud will lead to the formation of an accretion disk to 
conserve angular momentum; thus much, if not most, 
of the mass of low-mass stars is probably accreted from 
disks. However, if disks transport the infalling mass 
steadily to the central star, accretion luminosities are 
considerably higher than typically observed protostellar 
luminosities (Kenyon et al. 1990, 1994; Enoch et al. 2008; 
Evans et al. 2009). This "luminosity problem" can be 
alleviated if disks spend most of their evolution accret- 
ing slowly, undergoing episodic accretion outbursts where 
most of the mass is added to the central stars. One impli- 
cation of non-steady accretion is that mass may pile up 
in specific regions of disks, which could have significant 
consequences for understanding planet formation. 

Over the past decade, disk angular momentum trans- 
port by the magnetorotational instability (MRI) (Bal- 
bus & Hawley 1998) and by the gravitational instability 
(GI) (Durisen et al. 2007) has become much better un- 
derstood. It seems likely that both types of instability 
need to be considered to understand protostellar accre- 
tion. At early evolutionary stages the disk is likely to 
be quite massive with respect to the central protostar, 
suggesting that GI may be important; in addition, these 
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disks are so cold that thermal ionization cannot sustain 
the MRI, which requires coupling of the magnetic field 
to the mostly-neutral gas through collisions (e.g., Reyes- 
Ruiz & Stepinski 1995). On the other hand, external 
ionizing sources (cosmic and/or X rays) can provide the 
necessary ionization for MRI activitation up to a lim- 
iting surface density, ionizing the outer surfaces of the 
disk and leading to accretion in an "active layer" , leav- 
ing a "dead zone" in the midplane (Gammie 1996). In 
addition, at high accretion rates, the inner disk can be- 
come thermally-ionized, as in FU Ori outbursts (Zhu et 
al. 2007, 2009b). 

Mismatches in the transport rate between the GI and 
MRI can lead to outbursts of accretion (Armitage et al. 
2001, Zhu et al. 2009 be). In Zhu et al. 2010 (Pa- 
per I), we constructed one-dimensional, time-dependent 
disk models to study the protostellar unsteady accretion. 
These models assume both thermally-activated and lay- 
ered MRI-driven accretion along with a local treatment 
of the GI. We further adopted steady mass addition at 
a specified outer disk radius to drive the system. We 
found that accretion in protostellar disks is unsteady over 
a wide range of parameters because of a mismatch be- 
tween GI transport in the outer disk and MRI transport 
in the inner disk. Our results for outburst behavior in 
these one-dimensional model were sufficiently compara- 
ble to our previous two-dimensional simulations of out- 
bursts (Zhu et al. 2009c), confirming the utility of the 
faster ID simulations to explore wider ranges of param- 
eter space. 

In this paper we extend our results from Paper I 
to study long-term disk evolution from the protostellar 
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phase to the T Tauri phase, using a more self-consistent 
treatment for infall from a collapsing, rotating cloud core. 
This allows us to study the effect of the initial core rota- 
tion and different disk accretion configurations (layered 
accretion, Gl-only accretion, and constant-a accretion) 
on the disk formation and evolution. Rice, Mayo, & Ar- 
mitage (2009) have also investigated somewhat similar 
one-dimensional models but do not assume layered ac- 
cretion. 

Our models treat the same phase of evolution as 
Vorobyov & Basu (2005, 2006, 2007) investigated with 
two-dimensional models. They also argue that proto- 
stellar accretion is generally non-steady, but for a differ- 
ent reason, specifically the accretion of clumps created 
by gravitational instabilities. Kratter et al. (2009) have 
also investigated the protostellar accretion phase with 
three-dimensional numerical simulations. Although our 
treatment of the core collapse and the GI is much more 
schematic than used in the above investigations, we are 
able to treat radiative cooling more realistically, and can 
follow disk evolution to much smaller radii, where the 
(thermally- activated) MRI can become important. Fur- 
thermore, our two-zone disk model allows us to study 
the effects of different disk accretion configurations effi- 
ciently. The variety of disk structures predicted by our 
simulations with different initial cloud core rotations and 
disk accretion mechanisms can be tested by future EVLA 
and ALMA observations, which will help us to better un- 
derstand disk accretion and planet formation. 

We describe our one-dimensional, two-zone model with 
self-consistent infall in §2. In §3 we explore the behavior 
of the disk models with infall. In §4 we discuss the im- 
plications of our results, and summarize our conclusions 
in §5. We defer all derivations to the Appendix. 

2. 1D2Z MODELS WITH INFALL 

The one-dimensional two zone (1D2Z) disk model has 
been introduced in Paper I. It consists of two vertically 
averaged layers (the surface layer and the "dead zone") 
evolving independently. The surface density evolves ac- 
cording to the mass and angular momentum conservation 
equations in cylindrical coordinates, 



R 



(1) 



2TrR-(Z t R 2 ny — (M l R 2 n) = 2k—(R 2 W R4> , 1 )+2kK 1 

(2) 

where £ is the total surface density, Q is the angular fre- 
quency, Mi = —2nY>iRvi is the radial mass flux in the 
disk, Wftcf,^ — RYiii>id£l / dR, Vi is the viscosity, and i de- 
notes 'a' (active layer) or 'd' (dead zone). 2irgi(R,t) and 
The terms 2irAi are the mass and angular momentum 
flux brought by the infalling matter into the disk (Cassen 
& Moosman 1981). Assuming instantaneous centrifugal 
balance, 4 equations 1 and 2 can be simplified to 

*-^,^, + ^ 

4 The validity of this assumption with infall has been shown in 
Appendix B in Cassen & Moosman (1981). Instantaneous centrifu- 
gal balance is a good approximation in our simulations because the 
change of the central star mass on an orbital time is very small. 
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where M r is the sum of the mass of the central star (M c ) 
and the disk mass within R, using the approximation 
of the gravitational potential adopted by Lin & Pringlc 
(1990). The first term in equation 3 represents disk ac- 
cretion due to viscosity; the second term is due to the 
central star mass changing with time; and the third term 
is due to the infalling matter (Cassen & Moosman 1981). 
Since the infalling matter only falls onto the surface of 
the disk, we assume gd and A^ are 0, thus g a and A a 
can be written as g and A for short. The effect of the 
infalling matter onto the disk is limited to just two free 
functions: g(R,t) and A(R, t). 

The surface densities E a and exchange mass as de- 
scribed in Paper I to maintain E a < S^, where Y,a is 
the maximum MRI active layer surface density and is 
assumed to be a constant throughout the disk (Gammie 
1996). We solve equation 3 and 1 sequentially: Mi is cal- 
culated with the current disk temperature and £j, and 
then it is inserted into equation 1 to update £j at the 
next timestep. 

The disk temperature is determined by the balance be- 
tween the heating and radiative cooling. Here the exter- 
nal temperature T ext , which represents the heating effect 
of the irradiation from the central star is assumed to be 



rpi 

ext 



fL 



4irR 2 c 



(4) 



where L is the total luminosity of the star and changes 
with the central star mass by 



8M. 1 
9M + 9 



(5) 



which is a fit to the T- Tauri star birthline from Kenyon 
& Hartmann(1995). The change of luminosity has lit- 
tle effect on the disk evolution in our simulations, a is 
the Stefan-Boltzmann constant, and / is the coefficient, 
which accounts for the non-normal irradiation from the 
central star at the disk's surface. We assume /=0.1 in 
this paper. The viscous heating rate of the disk is 



Qvisc — 2 Wr^Q , 



(6) 



where the stress Wr^ — 3/2Xji/jf2 assuming Keplerian 
rotation. Both gravitational and magneto-rotational in- 
stabilities are considered. The anomalous viscosity Vi is 



(7) 



where c Si are the sound speed of the active layer and the 
dead zone, at = uq + olm and 



a Q 



-Q 



(8) 



The Toomre parameter Q is calculated using the disk 
central temperature (T d ) and the total surface density 
(S a +Sd). The form of aQ is motivated by a desire to 
make gravitational torques become important only when 
Q < 1.4, as indicated by global three-dimensional sim- 
ulations (e.g., Boley et al 2006). The parameter olm is 
set to be 0.01 when the disk temperature is above the 
preset critical MRI temperature for thermal activitation 
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(Tm=1500 K is assumed in this paper) or at the sur- 
face active layer. Since the rise time of MRI is about 
the orbital timescale (Stone et al. 1996), we assume that 
the MRI viscosity grows on the orbital timescale up to 
a M whenever the temperature is above the critical tem- 
perature (Tm). The effects of different values of 
Tm, and S a on outbursts are discussed in Paper I. The 
treatment of radiative cooling is the same as in Paper I. 

Recent observations suggest that the density structure 
of protostellar molecular cores, when circularly-averaged, 
can be approximated by static Bonnor-Ebert spheres 
(Alves et al. 2001; di Francesco et al. 2007; Ward- 
Thompson et al. 2007). Though real clouds have more 
complex structure (e.g., Benson & Myers 1989), result- 
ing in more complicated patterns of infall (Kratter et 
al. 2009), we adopt this simplified structure to limit the 
number of model parameters. This leads us to use an ap- 
proximate model of infall corresponding to that expected 
from a critical Bonnor-Ebert sphere, which exhibits sim- 
ilarities to that of the collapse of the singular isothermal 
sphere (Terebey, Shu, & Cassen 1984). 

Analytical and numerical studies ( Foster & Cheva- 
lier 1992; Henriksen et al. 1997; Gong & Ostriker 2009) 
have shown that the collapse of the critical Bonnor-Ebert 
sphere can be divided into two stages. In the first stage, 
the roughly constant-density inner region collapses to the 
central regions nearly simultaneously; then in the second 
stage, infall from the remaining cloud core (p oc r~ 2 re- 
gion 5 ) is at a rate similar to c 3 c /G (c sc is the sound speed 
at the cloud core temperature) of the singular isotheral 
sphere (SIS) collapse model (Shu 1977). We assume that 
our simulation starts at the end of the first stage of very 
rapid collapse of low-angular momentum material, so the 
central star initially has the mass of the inner flat re- 
gion of the Bonnor-Ebert sphere, and the collapse of the 
remaining cloud follows the singular isothermal sphere 
solution (Shu 1977). 

We use a two-component density profile to represent 
the Bonnor-Ebert sphere as in Henriksen et al. (1997). 
The two density components are 



p = p c at £ < V2 , 
2p c C 2 at V2 < £ < 6.5, 



(9) 
(10) 



where p c is the central density, r is in the normalized 
unit £ = r/(Cg C /47rG/9 c ) 1 / 2 , c sc is the sound speed of the 
cloud and p = 2p c £~ 2 is the density profile of the SIS 
model. £=6.5 corresponds to the critical Bonnor-Ebert 
sphere radius and the radius of the fiat density region 
£ = V2 corresponds to r lc =(c 2 sc /2irG pc) 1 ^ 2 ■ The total 
mass of each component is 



M(r < ric ) = 



M(r ic < r < r cr ) 



0.27cL 



G 3 / 2 p. 



1/2 



2.88c 3 



G 3 / 2 p, 



1/2 



(11) 



(12) 



The total mass with this simplified density profile is 
2>.lc 3 c /G 3 ^ 2 pV 2 , which is close to the mass of a real crit- 

5 In this paper, R denotes the cylindrical radius in the disk while 
r denotes the radius in the cloud core. 



ical Bonnor-Ebert sphere 4.2c 3 c /G 3 / 2 pc • Since the flat 
inner region has 10% mass of the r~ 2 region, we assume 
the flat region has 0.1 M and the rest of the Bonnor- 
Ebert sphere has 1 M . 

The rest of the 1 M cloud core beyond ri c collapses 
in a manner similar to that of the singular isothermal 
sphere, in which case the inner regions collapse first 
due to their shortest free-fall times and this free-fall 
zone extends outwards linearly with time (moC sc t, where 
m =0.975.) (Shu 1977; Terebey, Shu, & Cassen 1984). 

Thus, the mass infall rate is Mi„=moC 3 c /G, which is con- 
stant during the collapse. At time t, all the material col- 
lapsing to the center was originally at the cloud radius 
ro=(m /2)c sc t-|-r ic , where the addition of r ic (the radius 
of the flat density region in the Bonnor-Ebert sphere) is 
due to the fact that we define t=0 as the time when the 
flat region in the Bonnor-Ebert sphere has already col- 
lapsed to the central star. If the protostellar cloud core 
is initially in uniform rotation with angular velocity f2 c , 
the material falling in from different directions will have 
different angular momenta and arrive at the midplane at 
differing radii within the so called "centrifugal radius" 
(Cassen & Moosman 1981) 



R c = r. 



n 2 



GM C 



(13) 



where M c is the central star mass. Assuming M c in- 
creases with t , R c oc t 3 , similar to that found by Terebey, 
Shu, & Cassen (1984) for the singular isothermal sphere 
collapse. 

We use the Cassen & Moosman (1981) solution for the 
infall to the disk at a given time from a spherical cloud 
initially in uniform rotation. We start the infall from 
the transition radius r ic at t = 0. Then the mass flux 
g(R, t) and the angular momentum A brought by the in- 
falling matter for a uniformly-rotating singular isother- 
mal sphere is, 



g(R,t) 



Min 

4irR c 



A(R,t) = g(R,t)R 



GM C 
Rc 



Rc 

1/2 



-1/2 



(14) 



if R < R c , (15) 



and 



g(R,t) = 0, 
A(R,t) = if R > R c . 



(16) 
(17) 



where Mj„ = (moC 3 c /G) is the infall rate and R c is the 
centrifugal radius defined in equation 13. Inserting equa- 
tions 14-17 into equations 1-3, we can simulate the disk 
evolution under the infall. 

3. LAYERED ACCRETION WITH INFALL 

In all cases we start with M c = 0.1 M Q and a cloud 
mass of 1 M Q . We adopt a temperature of 20K for 
the envelope. This parameter enters in setting the infall 
rate M in =c s c 3 /G = 4.2 x lO^Mgyr" 1 ; thus the infall 
lasts 2.4 xlO 5 yr, roughly consistent with observations 
(Kenyon et al 1990, 1994; Evans 2009). This leaves only 
one parameter to vary, the initial angular velocity of the 
cloud Q c , which affects the maximum centrifugal radius 
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onto the resulting disk. For a singular isothermal sphere 
(SIS) core, 

where we have assumed that the central star mass at the 
end of infall is approximately the cloud core mass (M co ); 
thus 

2 



Rc 



12 AU 



lO-^rads- 1 



(19) 



Since K C max directly constrains the disk size after the 
infall, it is the most important parameter determining 
the initial disk properties. 

3.1. Fiducial model 

We take as our fiducial model one with £^=100 g 
cm~ 2 , ct!M=0.01, and fi c =10 _14 rads _1 so that R cm ax — 
12 AU. 

Figure 1 shows the fiducial model's mass accretion rate 
and the central star and disk masses as a function of time. 
The overall evolution can be divided into three distinct 
stages. Initially there is a phase of quasi-steady disk ac- 
cretion, with the accretion rate matching the infall rate 
(M in = 4.2 x 10 _6 M Q yr _1 ). This occurs because the 
infall is to small radii where the disk can become hot 
enough to sustain the MRI thermally. As infall proceeds, 
thermal activation of the MRI cannot be sustained, and 
gravitational instability takes over as the main mecha- 
nism of mass transfer. However, the GI is inefficient at 
these radii and so matter piles up. Eventually, the disk 
becomes opaque enough that the radiative trapping of 
the energy dissipated by the GI raises the central disk 
temperature sufficiently to thermally-activate the MRI, 
and an outburst of accretion occurs. As shown in Paper 
I, the result is cycles of outbursts transferring material 
inward; the GI moves material in slowly until enough 
mass is accumulated with enough viscous dissipation to 
trigger the thermal MRI outbursts. Finally, after infall 
stops, the disk is not massive enough to sustain GI, and it 
accretes only from the MRI active layer. Since the outer 
disk has a surface density smaller than the maximum 
surface density that can be ionized by cosmic and/or X 
rays, it becomes fully viscous due to the MRI; this leads 
to continued expansion of the outer disk, with a slow 
draining of material into the inner disk and central star. 
During the T Tauri phase, the layered accretion itself 
may also trigger outbursts, depending on the detailed 
layered structure (see below). 

The disk structure during the outburst stage is illus- 
trated in Figure 2, shown at the time when the centrifu- 
gal radius (R c ) is 10 AU (labeled by the triangle). Al- 
though the infall ends at R c , the disk extends far beyond 
R c due to the outward mass transfer by the active layer. 
With the continuous infall, the disk is marginally grav- 
itationally stable within R c , and the energy dissipation 
due to the GI has significantly heated the disk so that 
the outburst will be triggered at 2 AU (the upper right 
panel shows that T c is approaching 1500 K). At the time 
when infall ends, the region between ~1 AU and ^20 
AU is still marginally gravitationally unstable, although 
the disk within 1 AU has been depleted by the previous 
outburst (Figure 3). 



The innermost disk (R<0.1 AU) is purely MRI ac- 
tive due to the high temperature there produced by vis- 
cous heating and irradiation. Due to the efficient mass 
transport by the MRI, this inner region can be depleted 
rapidly, but is limited by the accretion from the outer 
active layer. Eventually, a balance is achieved, and the 
inner disk behaves like a constant-a disk with mass ac- 
cretion rate equal to the mass accretion rate of the in- 
nermost active layer (e.g., Gammie 1996). 

In the standard layered accretion model, the disk mass 
accretion rate is set as the active layer accretion rate at 
the dead zone inner radius where the MRI becomes ther- 
mally activated (e.g., Gammie 1996). We use the modifi- 
cation of the accretion rate including irradiation derived 
by Hartmann et al. (2006; see also the Appendix) 



M = 6.9 x 10" 



, R \ ( M c 
0.2AU J \IM Q 

1/4 



-1/2 



/ CUM \ 

uo-v 



lOOgcm 2 



T \ 1/4 / f \ 1/4 



where / is defined in equation 4. At the sublimation 
radius, f may be significantly larger than 0.1. With our 
numerical inner boundary 0.2 AU, the simulation shows 
the disk has a mass accretion rate 1O~ 8 M yr^ 1 in the 
quiescent state, which agrees with the above estimate. 

As discussed in Paper I, we find unstable behavior at 
the dead zone inner radius because of the rapid change 
in cxm ■ This accounts for the modest but significant vari- 
ability seen during the quasi-steady state 6 . This insta- 
bility cannot be reliably treated in our model, due to the 
complexity of the processes involved including MRI ac- 
tivation and dust sublimation; thus the precise nature of 
the variations in accretion during this phase is uncertain. 

A massive, marginally gravitationally unstable dead 
zone remains well after the end of infall (upper left and 
lower right panels in Figure 3). This dead zone is main- 
tained because the active layer mass accretion rate in- 
creases with radius (Gammie 1996). From the mass con- 
servation equation one can derive the time dependence 
of the total surface density E=£ J 4+E<j, 



as 

dt 



d 



RdR 




1/4 



J_ 

0.1 



M 

1/4 



-1/2 



10- 2 



lOOg cm 



g cm yr 



(21) 



Thus 9S /dt > and cx The disk surface den- 

sity thus increases linearly with time, more rapidly at 
smaller radii. When the surface density increases to the 
extent that the disk becomes gravitationally unstable, 
self-gravity transfers the excess mass and self-regulates 
the disk to an approximate Q~l state. From equation 
21 one can see that it only takes 10 5 yr for the disk at 10 

6 Although wc avoided this instability at the T Tauri phase by 
choosing the inner numerical boundary at 0.2 AU, this instability 
is present in our simulation in the first 10 5 years (short variations 
in Figure 1) when the dead zone inner radius is larger than 0.2 AU. 
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AU to become gravitationally unstable by layered accre- 
tion. Thus, active-layer accretion can maintain a dead 
zone that is marginally gravitationally unstable well after 
infall stops. 

A sharp density drop appears at the outer radius of 
the dead zone Rd (also shown in Reyes-Ruiz(2007) and 
Matsumura et al.(2009)) and Rd starts at the maximum 
centrifugal radius after the infall and gradually moves in- 
ward with time, so that the dead zone becomes increas- 
ingly narrow as the system evolves (see the Appendix). 
In the outer disk beyond Rd where only the active layer is 
present, the disk surface density evolves from a E a R^ 1 
distribution to an asymptotic E oc R distribution at 
large times. This result can be understood by modeling 
it as a pure a disk but with a fixed E inner boundary, as 
discussed in the Appendix. 

3.2. Varying core rotation 

The initial disk size is constrained by the maximum 
infall radius R cma x- We next consider the results for 
n c = 3 x l(T 15 rad s" 1 and 2 x KT 14 rad s" 1 . 

Figure 4 shows the accretion behavior of these models 
during the first 2 Myr. In the slowly-rotating case, all of 
the mass is accreted in a quasi-steady fashion, and the 
remaining disk mass is very small, whereas the collapse 
of the rapidly-rotating core results in very little quasi- 
steady accretion, an extended period of outbursts lasting 
well beyond the end of infall, and a very massive disk. 

The differing behavior of the three models can be un- 
derstood with reference to Figure 5, which is a modifica- 
tion of the M vs. R plane discussed by Zhu et al. (2009b). 
For the slowest-rotating core model, K cmax ~ 1 AU 
(equation 19). Reference to Figure 5 shows that this 
radius is just within the maximum radius Rm for steady 
MRI-driven accretion at M = 4 x 10~ 6 M Q yr" 1 . In this 
region, the disk can maintain a temperature greater than 
the critical value for MRI activation; the disk then ac- 
cretes on to the star as a typical a-disk, with modest 
oscillations due to the instability at the inner disk edge 
described in the previous section (see also Paper I). In 
the case of our fiducial model, R cm(a <~ 12 AU; this 
places the disk well beyond the region where the mis- 
match between GI and MRI occurs, resulting in out- 
bursts. The fiducial model accretes quasi-steadily for 
the first 8 x 10 4 yr; at this point R c ~ 1 AU. For the 
rapidly-rotating case, H cmax <~ 50 AU; very little mass 
infalls to radii within the MRI-stcady region, and most 
of the mass then gets accreted during outbursts. Because 
the mass transport is dominated by the GI at large radii, 
a relatively massive disk results. Note that the outer re- 
gions of the disk reside near the limit where gravitational 
fragmentation rather than simple transport might occur 
(Appendix; §4.4). 

Thus the range of rotations chosen here (given our as- 
sumed model for the initial cloud core) approximately 
span the extremes of disk evolution; at smaller ftc, the 
accretion is almost spherical and occurs quasi-steadily, 
while for larger values of f2 c gravitational fragmentation 
is likely to occur, considerably modifying disk evolution. 

The disk mass and surface density distribution are also 
significantly affected by the initial core rotation. With 
O c = 2 x 10~ 14 rad s _1 (left panel in Figure 6), although 
the dead zone shrinks from 35 AU at the end of the infall 



phase (0.3 Myr) to 10 AU at the T Tauri phase (1 Myr), 
the disk is very massive (0.3 M©) and still gravitationally 
unstable within 10 AU during the T Tauri phase. How- 
ever, for smaller rotation fl c = 3 x 10 _15 rad s _1 (right 
panel), the dead zone is significantly smaller (2 AU) and 
has less mass (0.05 M©). 

When the initial core rotates slowly and R C mox is small, 
the dead zone size is determined by R CTOO x- However, if 
the inital core rotates fast and R cm ax is large, the dead 
zone size is constrained by the radius where the active 
layer becomes gravitationally unstable, as discussed fur- 
ther in the Appendix. 

One important point is that even though the dead zone 
is small < 50AU, it can contain a significant fraction (or 
most) of the disk mass. Figure 7 shows that, no matter 
what the infall it is, most of the disk mass can be at small 
radii, rather than the large radii probed by current mm 
and sub-mm observations (in comparison with constant- 
ly models). 

3.3. Different disk configurations 

The ionization structure of disk surface layers depends 
on the dust size spectrum and the flux of ionizing ra- 
diation and cosmic rays, and is therefore poorly con- 
strained. Angular momentum transport in magnetically 
decoupled disks is also poorly understood, although re- 
cent work by Lesur & Papaloizou(2009), following Pe- 
tersen et al.(2007), strongly suggests that a nonlinear 
instability driven by baroclinicity and radiative diffusion 
may give rise to hydrodynamic turbulence and angular 
momentum transport. Therefore we have investigated 
disks with a reduced active layer (§3.3.1), and a residual 
dead zone viscosity (§3.3.2). 

3.3.1. Disks with less active layer and pure GI disks 

Little is known about the disk active layer, which is de- 
pendent not only upon cosmic and/or X ray fluxes, but 
on the presence of dust which can absorb charges and re- 
duce the MRI or even prevent it (Sano et al. 2000). Am- 
bipolar diffusion at the low density active layer beyond 
10 AU may further depress the effectiveness of the MRI 
(Chiang & Murray-Clay 2007). In addition, there may be 
no dynamo in the outer disk. Thus we consider disks with 
lower active-layer surface densities (E^ < lOOg cm -2 ). 

The long-term disk evolution is determined by the 
product of oim^a (equations 20 and 21); thus, as long 
as the product Y>aolm stays the same, disk evolution is 
similar. Figures 2 and 3 show disks with E^ = 100 g 
cm~ 2 and qm = 0.01 and disks with E^ = 10 g cm~ 2 
and ctu =0.1 have similar surface density, temperature, 
and mass accretion rate distributions at all times, only 
except that a lower E^ leads to a less massive outer disk. 

If um^a is lower, the disk accretion rate is lower than 
our fiducial model; thus the infall mass takes longer to 
accrete onto the star. In addition, with a smaller E^, the 
dead zone outer radius (Rd) is larger (equation A13). 

In the extreme case when there is no active layer 7 , the 
disk evolves towards Q~l during the simulation (Fig- 
ure 8), since the GI is the only mechanism to transport 
mass from the outer disk inwards. The disk is very mas- 
sive and has a sharp outer disk edge roughly at R cm0 z; 

7 We allow thermal activation of the MRI when the disk tem- 
perature is higher than Tjk- 
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(Vorobyov & Basu 2005, 2006, 2007). Without the ac- 
tive layer replenishing the innermost disk and the inner 
disk becoming gravitationally stable (due to the strong 
irradiation from the central star) , the disk accretion rate 
drops to zero (upper left panel in Figure 8). During the 
T Tauri phase, the disk evolves to Q~2 and stays at this 
gravitationally stable state (Figure 8). Thus the disk 
is massive (0.2 M Q ) during the T Tauri phase. Even if 
the disk has some small residual viscosity at Q>>1, it 
has limited effect on disk evolution simply because, with 
a < 10~ 4 , the viscous timescale is too long (more than 2 
Myr to transport mass from 50 AU to the central star) 
(see also Vorobyov & Basu 2009). 

3.3.2. Residual dead zone viscosity 

In our fiducial model, we assume the dead zone has 
no viscosity at all. However, numerical MHD simula- 
tions have shown that the dead zone may not be totally 
"dead", with a small residual viscosity. The viscosity 
could be due to the operation of the MRI with the ions 
turbulent mixing from the active layer (Ilgner & Nel- 
son 2008) or from turbulence penetrating the dead zone 
from the active layer (Turner et al. 2007). Recent work 
by Lesur & Papaloizou (2009) also persuasively demon- 
strates the existence of a nonlinear instability driven by 
baroclinicity and radiative transport (the magnitude is 
not yet clear, and the existing evidence suggests that the 
instability weakens when radiative diffusion is small — 
when the Peclet number is large — as might be expected 
in a dead zone with large optical depth. 

If the residual viscosity in the dead zone is a r d ~ 10~ 5 ), 
it has a negligible effect on disk evolution. Disk structure 
evolves simply on the viscous timescale; for a = 10~ 5 , 
this is long compared to T Tauri lifetimes ~ few xlO 6 
yr- 

However, if the residual a. r d >10 , the dead zone 
can transport the mass at a rate comparable to the ac- 
cretion rate of the active layer, with Mj cx £,a, and 
Yid » S a . We assume the disk has a constant mass 
accretion rate after a long period of time, and at the 
stellar surface dfl/dR=0. Furthermore, if the disk is 
irradiation-dominated such that TocR -1 / 2 , ^ a =k a R and 
i y d=k ( iR where k and k^ are constants and are propor- 
tional to a a and ad, we derive (in the Appendix) 



M 



3-n-kdR 



1 - 



R* 



1/2- 



— — E 

kd 



(22) 



If the dead zone transports mass at a rate higher than the 
active layer, the last term can be neglected for an order of 
magnitude argument. Compared with a constant-a disk 
having the same mass accretion rate, the layered disk 
has a surface density higher than the constant-a disk by 
a factor of a/ ad- If we assume the outer pure MRI ac- 
tive disk and the inner layered disk has the same mass 
accretion rate, and the outer disk only has v a while the 
inner dead zone is dominated by Ud, the disk's surface 
density will change by a factor of a a /ad within the dead 
zone outer radius Rd, which is shown in Figure 9. Also 
this density change at Rd is more gradual (equation 22 
is a smooth function of R) than the zero dead zone vis- 
cosity case. This can have a significant impact on planet 
migration (§4.6). 



4. DISCUSSION 
4.1. Dead zones and star formation by disk accretion 

The main result of our work is that star formation 
through disk accretion, assuming a sufficiently large ini- 
tial disk that the MRI is not thermally activated through- 
out the disk, generally results in a dense belt of material 
- the "dead zone" -which has implications for disk struc- 
ture and planet formation. 

To demonstrate this more clearly, we show the results 
for a constant-a disk for the same infall model. By 
comparison with layered disks in Figure 7, constant-a 
disks have similar outer disk but lack the massive dead 
zone. Constant-a disks respond to infall just as layered 
disks do: a faster rotating core leads to a more massive 
disk since the infall mass onto larger radii takes longer 
to accrete to the central star. Our models show that 
with the same a in constant-a models, cores with O c = 
2 x 10~ 14 rad s _1 lead to disks 10 times more massive 
than the disks produced by cores with fl c = 3 x 10~ 15 rad 
s _1 at 1 Myr. However, even if we allow constant-alpha 
disks to have different values of alpha (Figure 10), all 
these models produce surface densities far in excess of 
100 g cm" 2 over large areas of the disk at the end of 
infall (lower left panel of Figure 10). Thus, a reduced 
angular momentum transport efficiency associated with 
the failure of the MRI at S > 100g cm— 2 implies the 
formation of a dense belt or dead zone in the disk. 

Another feature of the constant-a model is that, al- 
though the core rotation determines the disk mass, the 
disk surface density shape depends most strongly on a. 
This can be simply explained by the similarity solution 
of the constant-a disk where the shape of the disk is de- 
termined by the scaling radius (Hartmann et al. 1998), 
which is proportional to t and a. Thus as long as t and a 
are the same, the disk surface density has similar shape. 
We also notice that, even for a constant-a disk, the disk 
is gravitationally unstable at the infall phase with a rea- 
sonable set of parameters (a=0.01 and O c =2x 10~ 14 rad 
s" 1 ) (Figure 10). 

4.2. Relation to previous work 

As mentioned in the introduction, previous investiga- 
tions by Vorobyov & Basu (2006, 2007) and Kratter et 
al. (2008, 2009) have considered the evolution of disks 
formed by infall from rotating cloud cores with angular 
momentum transport by GI. Vorobyov & Basu find out- 
bursting behavior, as we do, but for different reasons; 
specifically, in their models time-dependent accretion is 
driven by gravitational clumping, which we cannot treat 
in one dimension. However, with more realistic energy 
equations, it is harder to lead to strong gravitational 
clumping. Kratter et al. (2009) studied more generalized 
models of infall; they did not find clumping but found 
conditions for which fragmentation would occur, which 
agrees with our estimate in Appendix B. 

Unlike these investigations, we study slowly rotating 
cloud cores forming single stars without disk fragmenta- 
tion, focussing on disk evolution after the central core 
collapse. We treat the GI as a local phenomenon, which 
may be problematic if the disk is very massive (Lodato & 
Rice 2005). However, during later stages when the disk is 
less massive, a local treatment appears to be reasonable 
(Lodato & Rice 2004, Cossins et al. 2009a). In addition, 
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we are able to treat radiative cooling more realistically; 
the treatment of thermal physics plays a crucial role in 
GI transport (e.g., Durisen et al. 2007) and MRI-GI out- 
bursts (Paper I). 

Lodato & Rice (2005) found recurrent episodes of 
gravitational instability in a massive disk. Boley & 
Durisen (2008) considered three-dimensional simulations 
of gravitationally-unstable disks with accurate radiative 
transfer, and also suggested that high mass fluxes could 
result from the rapid onset of GI, producing something 
like an FU Ori outburst. The question is whether this 
behavior is a transient due to an initial condition, which 
may or may not be realizable as part of a natural evolu- 
tionary sequence. This GI outburst could be a boost in 
the MRI activation. 

Rice et al. (2009) also considered the evolution of disks 
formed by infall with GI as the transport mechanism. 
They also find that massive disks result, which leads 
them to suggest that another mechanism of viscosity 
must be operating to drain disks. We address this possi- 
bility in §4.5. 

Our treatment differs from all of the above in includ- 
ing MRI transport with layered accretion. This results in 
our outer disks tending to evolve viscously rather than 
through GI, so that mass accretion can continue even 
after the disk is gravitationally stable. Perhaps most im- 
portantly, our one-dimensional model allows us to treat 
the very innermost disk at the same time considering the 
outer disk evolution, something that is extremely difficult 
to do with two- or three-dimensional codes. This allows 
us to demonstrate behavior not seen in the other simu- 
lations; specifically, the onset of MRI-driven accretion in 
the inner disk, which tends to dominate the outburst be- 
havior - or, as we have shown, result in a phase of quasi- 
steady disk accretion while low angular momentum core 
material is infalling. 

At high accretion rates, it appears almost certain that 
the MRI should operate, as temperatures become high 
enough to sublimate dust, which might otherwise absorb 
charges. This notion is supported by observations of FU 
Ori (Zhu et al. 2007, 2008), which indicate that the disk 
can be hot enough to eliminate dust out to radial scales 
of order 1 AU. Thus we argue that inclusion of the MRI 
is important for any understanding of accretion onto cen- 
tral protostars. 

The behavior of non-thermally activated MRI by X- 
ray or cosmic ray radiation is much more uncertain, but 
its inclusion can also have important effects not seen in 
pure GI treatments. In our layered simulations, during 
the T Tauri phase the GI has limited effects on the long 
term disk evolution and the disk mass accretion rate is 
controlled by the active layer, as discussed in §3. To 
confirm this, we have run the same simulation as before 
but with the cxq from Armitage et al. (2001), which is 
based on Lin& Pringle (1987, 1990). We found that the 
different forms of olq have little effect on the long term 
evolution of a layered disk; the T Tauri disk evolution is 
more determined by the active layer. In both treatment 
of olq , the layered accretion leads to a relatively massive 
"inner belt" within 30 AU and a viscous disk beyond 30 
AU (Figure 3). 

4.3. Outbursts and luminosity problem 



Our model predicts that protostars will undergo FU 
Ori-like outbursts of rapid mass accretion to accumulate 
a significant amount of their mass, helping to solve the 
luminosity problem (§1). Here we look at the question 
of luminosity evolution more closely. 

The largest current survey of star-forming regions 
bearing upon this issue is that of the c2d Spitzer Legacy 
project (e.g, Enoch et al. 2009, Evans et al. 2009), 
which strongly suggests that protostellar accretion is 
time variable, with prolonged periods of low accretion 
rates. Evans et al. (2009) estimate that half of the mass 
of a typical low-mass star is accreted during only ~ 7% of 
the Class I lifetime, which they estimate as ~ 0.5 Myr. 
In our fiducial model, roughly half of the mass is ac- 
creted quasi-steadily during the first <~ 10 5 yr, with the 
other half being accreted during outbursts within the 
next <~ 1.5 x 10 5 yr. This agrees with the c2d obser- 
vations that half of the mass is accreted during high 
mass accretion rate stage. However, the outburst behav- 
ior (episodic accretion) for the individual protostar does 
depend on the initial rotation of the cloud core. The 
O c = 2 x 10~ 14 case accretes almost entirely via out- 
bursts, while the low- rotation case accretes essentially 
all of the stellar mass quasi-steadily. 

Evans et al. (2009) found that 59% of Class 0/1 stars 
have bolometric luminosities smaller than 1.6 Lq. How- 
ever, in our fiducial model, 30% of the Class 0/1 phase 
is in the quasi-steady accretion phase with luminosity 
-13 L (GM^M/R* by assuming M=4xlO" 6 M© y^ 1 
and M*/R*=M /1OR Q ), much larger than 1.6 L Q . One 
possible solution of this problem is that the extinction 
toward the protostar during this quasi-steady phase is 
high, with Ay > 400 using our model parameters; it 
might thus be difficult to identify these objects in such 
early stages of evolution. 

Another possible solution is the oscillation of the ther- 
mal MRI region (Wiinsch et al. 2006) can lead to addi- 
tional variations in the mass accretion rate and accretion 
luminosity (which is the reason we call the first steady ac- 
cretion phase "quasi-steady" ; variations are clearly seen 
in Figure 1 during the first 10 5 years) (see §3.1). We 
have avoided concentrating on this feature because it de- 
pends on the complex physics at the inner edge of the 
dead zone, which is unclear to us now. Nevertheless, 
this potentially unstable behavior may be important in 
understanding protostellar accretion in general and may 
even provide a mechanism for explaining EXor outbursts 
(Herbig 1977). Further investigation of this problem with 
a realistic calculation of the onset of the MRI is needed. 

Some other limitations of our models are that the in- 
fall is axisymmetric with a constant infall rate, and that 
the GI is treated as a local effect. If the infall has a 
complex pattern and the GI has a strongly global char- 
acter, the angular momentum transport during the infall 
phase could be very efficient (Larson 1984) and disk ac- 
cretion may be more variable (Vorobyov 2009) or even 
form warped disk (Bate et al. 2009). 

There is also a distinction between the model pre- 
dictions for outbursts during the infall phase and those 
which occur after infall ends, during the T Tauri phase. 
Though they are both triggered by mass accumulation, 
this is caused by GI during the infall phase, which is in- 
dependent of the layered structure (Paper I), while the 
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outbursts during the T Tauri phase are triggered by the 
mass accumulation from active layer accretion, which is 
sensitive to the assumed value of S . 

4.4. Rotation and Fragmentation 

Our most rapidly-rotating model formally predicts FU 
Ori-like outbursts during the T Tauri phase, driven by 
active layer accretion. However, this T Tauri phase out- 
burst behavior is uncertain because it sensitively de- 
pends on active layer properties and the massive disks 
which would produce such behavior fragment instead into 
multiple stellar systems with individually smaller disks 
(Rapidly-rotating cores tend to produce large disks that 
are subject to gravitational fragmentation at R> 50 AU 
in a realistic disk, see Rafikov 2009; Cossins et al. 2009b; 
see Appendix). 

Finally, we consider possible values of envelope angu- 
lar momentum. Typical observational estimates of the 
ratio of rotational kinetic to gravitational energy span 
a range of three orders of magnitude (e.g., Caselli et 
al. 2002), and thus are difficult to apply to our simu- 
lations. Moreover, the total angular momentum is ex- 
tremely sensitive to assumptions about the distribution 
of mass with radius as well as the (not necessarily uni- 
form) angular velocity; as cores are generally not Bonnor- 
Ebert spheres, or even axisymmetric about the rotation 
axis (Tobin et al. 2009), making it difficult to compare 
with observations. As mentioned above, cores with more 
angular momentum than our fastest-rotating model are 
likely to fragment into multiple systems. Whether sig- 
nificant numbers of slowly-rotating cores exist which can 
produce diskless T Tauri stars, or T Tauri systems with 
very low-mass disks, will require more detailed obser- 
vational analyses of the complex geometries of starless 
cores. 

4.5. Dead zones: do they exist? 

We have shown that a natural consequence of our lay- 
ered model is the formation of dead zones that may 
remain relatively massive for more than 10 6 yr (also 
Vorobyov & Basu 2006, 2007; Rice et al. 2009). This 
stands in contrast with purely viscous models of disk 
evolution, as shown in Figure 7. An advantage of highly- 
viscous models is that they provide an explanation for 
disk accretion rates during the T Tauri phase (Hartmann 
et al. 1998; Dullemond, Natta, & Testi 2006), whereas 
the layered accretion model has some difficulties in ex- 
plaining T Tauri accretion rates (Hartmann et al. 2006). 
On the other hand, the pure viscous models are unable 
to explain FU Ori outbursts, which require large mass 
reservoirs at distance scales on the order of an AU (Ar- 
mitage, Livio, & Pringle 2001; Zhu et al. 2007, 2008, 
2009). 

At present the best that can be said is that the MRI 
is likely to operate at r < 0.1AU, and it is unlikely to 
operate at larger radii where the temperature is lower. 
Even the standard form of the minimum-mass solar neb- 
ula (Weidenschilling 1977), let alone the recent amplifi- 
cation suggested by Desch (2007; Figure 12), indicates 
surface densities well in excess of those thought to be 
penetrated by cosmic rays, or even stellar X-rays. At any 
rate, the associated change (presumably decrement) in 
angular momentum transport efficiency where the MRI 



shuts off is likely to give rise to a feature in the surface 
density. 

At first sight, the surface densities of the dead zones 
in our evolutionary models seem unreasonably large in 
comparison with previous models. However, Terquem 
(2008) has shown that, even for a constant mass accre- 
tion rate layered model, the dead zone is very massive. 
Reyes-Ruiz (2007) has also shown that a massive disk 
still exists even if the dead zone has some residual viscos- 
ity. Some of our models do not predict values of £ much 
higher than that of the Desch (2007) model (Fig. 12). If 
the dead zone has some residual viscosity with a > 10~ 4 , 
the dead zone mass can be significantly reduced but still 
relatively massive compared with a constant-a disk (Fig- 
ure 13). Furthermore, some exoplanetary systems have 
masses substantially larger than that of Jupiter, suggest- 
ing that more than minimum mass solar nebulae are re- 
quired, especially if migratory loss of solid bodies is sig- 
nificant. 

Have high surface density regions been ruled out by 
recent interferometric observations of submm dust emis- 
sion in disks, such as those of Andrews et al. (2009)? To 
examine this possibility, in Figure 13 we compare various 
models (our fiducial model, fiducial model with a residual 
dead zone viscosity, and a constant-a accretion model) 
with the suite of disk structures inferred by Andrews et 
al. (2009). The agreement between the Andrews et al. 
results and our models in the outer disk is not surpris- 
ing, because they adopted viscous disk solutions like ours 
to model the observed interferometric visibilities. If we 
adopt the same opacity as Andrews et al. , we find that 
there is no strong signal in the disk surface brightness 
distribution at the edge of the dead zone because it is 
optically thick and the temperature distribution is domi- 
nated by irradiation from the central star (the brightness 
temperatures for the layered model and the constant-a 
model are the same at 1 mm in Figure 14). The model 
does predict a possible feature at A ~ 1 cm where the 
outer disk becomes optically thin (Figure 14); this might 
be testable with the EVLA. The high surface densities 
in the dead zone may lead to rapid coagulation of the 
solid particles, however, reducing the optical depth con- 
trast between the dead zone and the MRI-active outer 
disk region. 

Our results suggest the potential importance of dead 
zones for simply providing a large mass reservoir for solid 
body formation. Higher density regions may be prob- 
lematic in that timescales for radial drift of small bod- 
ies and Type I migration become much shorter; on the 
other hand, this may be offset by having more material 
to start with, and other features of the dead zone may 
help provide "traps" where planetesimal can form effi- 
ciently (Rice et al. 2004) and migration can be halted or 
even reversed (see the following section). 

In this connection it is worth noting the existence of 
the so-called "transition" disks, where outer, optically- 
thick disks surrounding T Tauri stars have very large 
inner disk holes - either partially or fully evacuated (e.g, 
Espaillat et al. 2007). In some cases these large inner disk 
holes may be the result of tidal torques from companion 
(binary) stars, but in others giant planets may be the 
cause of the inner disk clearing. In the latter case, it 
may be necessary to have multiple giant planets form 
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nearly simultaneously over a range of radii in order to 
explain the sizes of these disk holes. The dead zones of 
our models, with relatively large surface densities over 
a significant range of radii, may be able to promote the 
necessary runaway growth. 

4.6. Rd and planet migration 

The sharp density jump at the dead zone outer radius 
(Rd) can significantly affect the planet migration outside 
the dead zone (Matsumura & Pudritz 2007); the associ- 
ated torque on the planet depends on Ra, the density 
jump factor F, and the width of the density jump. 

In the framework of this paper, we derive the first two 
parameters analytically and test it numerically. For Rd, 
§3.1 and the Appendix have shown that it starts at the 
outer disk at the end of the infall, and then travels in- 
wards with time at a speed given in equation A16 in the 
Appendix. With the initial R_D^R C mair (if we om y con- 
sider the slow rotating core: no GI fragmentation and 
active layer GI as discussed in the Appendix) and the 
speed of Rd, the position of Rd can be derived at any 
time later (Figure 11). Based on equation A16, the speed 
of Rd gets slower when moving inwards, which is seen in 
Figure 11. 

We can derive F if the dead zone residual a is small 
enough (a < 10~ 5 ). With such a small dead zone resid- 
ual viscosity, the disk is gravitationally unstable within 
Rd- Thus the surface density jump can be simply de- 
rived by dividing the surface density of the Q=l disk to 

Sq=i = c s n = /45A[A 7/4 / 100 gC m- 2 \ 

E A ttGEa \ Rd ) V Z A ) 

(23) 

Here we have assumed the temperature profile as T=200 
K(R/1 AU) -1 / 2 . Thus, as the dead zone moves inwards 
(Rd becomes smaller), F increases. However, the den- 
sity jump width, which is the last parameter required to 
calculate the torque on the planet, cannot be constrained 
by our ID simulation. 

4.7. MMSN and planet formation in the layered disk 

Although at the early stage the inner disk is massive, 
at later stage (~Myr) the outer disk (beyond 10 AU) is 
comparable to the minimum mass solar nebulae (MMSN) 
from Weidenschilling (1977) (Figure 12). Due to the 
boundary effect at the dead zone outer radius Rd (§3.1 
and Appendix), the outer disk evolves towards E cx R 
as in the standard MMSN. Furthermore, if the dead zone 
is massive, Rd moves inwards very slowly (§4.5); then 
E cx R~ 15 lasts for a long time. 

If planets form in a massive dead zone, they may be lost 



by inward migration; however, some may be trapped at 
the inner boundary (Kretke et al. 2009) or outer bound- 
ary (Matsumura et al. 2007,2009) of the dead zone. 

5. CONCLUSIONS 

In this paper, we have constructed a one-dimensional 
two-zone accretion disk model to study disk formation 
and long-term evolution under the collapse of a BE ro- 
tating core. The model evolution can be divided into 
three stages. At the early stage, when the mass falls to 
the inner disk within AU scale, the MRI can be sustained 
in the inner disk and efficiently and steadily transfers 
the infalling mass to the central star. Later, when the 
mass falls beyond AU scale, the disk goes to the outburst 
stage due to the accretion rates' mismatch by the MRI 
and GI as described in Paper I. After the infall com- 
pletes, the disk enters the T Tauri phase and evolves on 
its own. Cores with higher initial rotation end up with 
a more massive disk and more disk episodic accretion 
events (outbursts). As long as the initial cloud core does 
not rotate extremely slowly to form a tiny disk (R C max^ 
AU) , more than half of the star mass is built up by out- 
bursts, which eases the "luminosity" problem. 

Disks exhibit a variety of behavior during the T Tauri 
phase. For a disk with accretion sustained only by grav- 
itational instability, the disk evolves towards a Q=l disk 
and the disk truncates at a radius slightly larger than the 
maximum centrifugal radius of the infall. If the disk has 
an active layer at the surface, however, the active layer 
can extend to a much larger radius and a sharp density 
drop develops at a characteristic radius Rd that sepa- 
rates the marginally gravitationally stable dead zone and 
the MRI active but gravitationally stable outer disk. The 
density jump at Rd may be observable by the EVLA and 
ALMA. The formation of a dense belt of material is asso- 
ciated with the failure of magnetically driven transport 
due to low ionization at intermediate radius in the disk; 
the only ways to avoid this are (1) if there is a separate, 
equally efficient hydrodynamic transport mechanism, or 
(2) if for some reason the MRI fails in the outer disk as 
well, perhaps due to dynamo failure. 
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APPENDIX 
A. LAYERED ACCRETION 

For a layered disk dominated by viscous heating, Gammie (1996) showed that the central temperature and the mass 
accretion rate at the radius R are 

/ y a \ 4/5 

x 5) K , (Al) 

V200gcm- 2 y ' v ; 
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9/5 

Moyr- 1 , (A2) 



where the Bell & Lin (1994) opacity has been used. 

In the irradiation dominated limit, the disk temperature is the external temperature (equation 4), ctT 4 = fL^/AirR 2 . 
In this case the disk temperature and mass accretion rate with radii are 



M = 6.9 x 1CT 8 



/ R \ ( M c \ 1/2 / a M \ 
VlALV \IMq) uo-v 

(20i^)(l|) 1,< (ol) ,/4M ^ 1 - (A4) 

Comparing equations A4 and A2, we see that the disk is irradiation dominated if the active layer has surface density 
(XU < 10 2 g cm -2 ), or qm is smaller than 10~ 2 , or the luminosity is significantly larger than the stellar radiation 
(L>10L Q ). 

A. 1 Dead zone 

In either of the above cases the layered disk accretion rates increase nearly linearly with radius, which results in 
piling up mass in the dead zone at small radii. Using the mass conservation equation 

r m = ~m [to) ' (A5) 

<9£ jdt oc R~ x . Thus the layered disk surface density increases linearly with time and it increases more rapidly at 
smaller radii. If we assume the dead zone has zero residual viscosity (non-zero residual viscosity has been discussed in 
§3.3.2) and active layer viscosity ^ a =k_R", disk evolution becomes: 

^ = 3n ( n + I) kZ A R n - 2 . (A6) 



x 



dt V 2 . 
Thus 

T,(R,t) = 3n (n+ fcE A i?"- 2 t + C, (A7) 

which increases linearly with time due to the layered accretion, as shown by Matsumura & Pudritz (2007). 

The dead zone starts at 0.1 AU and ends at the radius where the disk surface density is smaller than Ha (the dead 
zone outer radius R_d), leaving only the MRI active layer ionized by cosmic and/or X rays at larger radii. This pure 
MM active outer disk beyond Rd evolves like a constant-a viscous disk. However since the active layer of the inner 
disk accretes mass inwards, the dead zone is gradually depleted so that moves inwards with time. 

To determine the dead zone size (Rd) during the layered disk evolution, we need to know its initial position just 
after the infall. Initially Rd should be close to R cma x inside of which the infall mass lands. However this is only true 
if the initial core rotates slowly and R C max is small. If the inital core rotates rapidly and R cm ax is large, the dead zone 
size is constrained by the radius where the active layer becomes gravitationally unstable. Generally, the dead zone 
cannot extend to R>50AU if S J 4=100 g cm -2 due to the active layer GI. This is because the GI can be very efficient 
in transporting mass and angular momentum when Q < 1, and the surface density with Q=l can be considered as an 
upper limit that the disk surface density cannot exceed during the evolution. For a layered disk dominated by viscous 
heating, Gammie (1996) has shown that the Toomre Q parameter at R is 

At the outer dead zone radius (Rd), £ = Sa; with the condition that Q>1, 
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For an irradiation dominated disk, we derive a similar condition: 

Q = ^Tz7 } - < A10 > 

assuming 

T(i?) = 200K(^ J )- 1 /2. (All) 

Then, 

Since Q(R D )>1, 

RD<is (m^ELl)"\ v . (A13) 



In either case, Rd < 50 AU. To test this, we computed a case with a fi c = 2 x 10~ 13 rad s _1 core. The core mass falls 
into Rcmax as large as 400 AU, but the dead zone beyond 30 AU is quickly cleared by the GI so that Rd is smaller 
than 30 AU during the T Tauri phase. 

After we know the initial Rd position, the motion of Rd can be derived by considering the mass conservation at 
Rd- For the region extending from Rd-A_R to Rd, 

+ (A14) 

Assuming the mass flux (or vr) at Rd is zero, which will be justified in some cases later, and considering the time for 
the dead zone within A R to be depleted is At, we got 

ife-^^)=0, (A15) 
At AR K ' 

AR S a M(R D ) 

VD= AF =VR ^-r2,R(x a+ z d) - (A16) 



thus 



where M is as equation A4. If we insert M we can derive vd« l/(S a + Considering £ a +£d is larger at smaller 
radii, the speed of Rd decrease with time when it is moving inwards. Finally, when Rd moves to the very inner radii 
where S<j » S , Rd is almost halted. 



A. 2 Outer pure MRI active disk beyond Rd 
Beyond Rd the disk is purely MRI active with a constant olm- Its surface density evolution can be solved by 

as _ 1 d 
~dt ~ RdR 



3i?V2 IvVR 1 / 2 



(A17) 



dR 

with the boundary condition J^h d —'Ea- If v — kR n , this solution can be simplified by dividing it into two parts 

S(i2,t) = S'(i2,t) + SA(^y) ' , (A18) 

where the first term on the right comes from the disk evolution with a zero surface density boundary condition at Rd 
(£#£,=0), and the second term is the effect of the non-zero boundary surface density at Rd- The first term behaves 
like the similarity solution with a R _1 part and a power law decrease part at larger radius, and it eventually decreases 
to zero with time. The surface density evolution is determined by the competition between the first and the second 
term. During the infall phase, the disk is massive, and the first term is larger than the second term, thus the disk 
behaves like a similarity solution. Then after infall stops, as the disk accretes and Rd moves inwards, both of these two 
terms decrease, and as shown above, Rd moves inwards more slowly. Eventually Rd can be considered independent of 
t; the first term goes to zero and the second term dominates so that we have S oc R (if n=l), which distinguishes 
it from the self-similarity solution. 

Next, we consider the late stage when Rd is independent of time analytically by the Green's function method. If 
the viscosity v is a function of radius, equation A17 is a linear equation for E and can be solved by a Green's function 
(Lust 1952; Lynden-Bell & Pringle 1974). If v = kR n , the above equation becomes 
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and the radial dependence of £ is a linear combination of the Bessel functions J +fl and J_ M , where fj, 2 = 1/4(2 — n) 2 
(Lynden-Bell & Pringle 1974). However for a fixed £ inner boundary condition at Rd, 

Z(R = R D ,t)=Z A , (A20) 

we can substitute £ with £' 

£(#, t) = £'(#, i) + T, A {-j^J . (A21) 



so that equation A19 and the boundary condition A20 changes to 

as 7 _ i a 

and 



(A22) 



Y,'(R = R D ,t) =0, (A23) 

which becomes a normal disk evolution equation of £' with a zero surface density boundary condition. The solution 
of E' is well studied and this disk expands with limit t — > oo, E' — > 0, thus from equation A21, we derive t — > oo, 
E -4- Ea(-R/-Rz5) - ™ -1 ^ 2 . Thus as long as the disk evolves long enough the impact of any initial condition will be 
washed out, and the boundary term dominates the surface density distribution with E oc i? - " -1 / 2 . 

In the irradiation dominated case, n=l, we transform equation A22 by writing x = R 1 / 2 and a = E'i? 3 / 2 (Pringle 
1991) to get 

da 3kd 2 a , A „„, 

1 /2 

with the boundary condition a = 0atx = .T£> = . The general solution is then 

/oo 
Axe~ x2t sin[\{x - x D )/c]d\. (A25) 
-oo 

where c 2 = 3fc/4. A\ is determined by the initial conditions. Following Pringle (1991), in order to obtain the Green's 
function, we set the initial condition with 

a(x,t = 0) = a S(x-xi). (A26) 
With the delta function Fourier transform J exp{2irx\i)d\ = S(x) and x > xd, we derive 

A x = -—sin[\(x D -x 1 )/c\. (A27) 

7TC 



Thus the solution for the initial condition A26 is 

2tt 1 /2 c 



a(x,X!,t) = , {exp[-(x - xi) 2 /3kt] 



— exp[— (x + x\ — 2x D ) 2 /3kt}}, (A28) 
Finally, the solution a(x,t) for any initial condition a(x 1 t = 0) = <j'{x) is 

*(x,t) = J ^^i/o^ 72 {exp[-(x - Xl ) 2 /3kt] 

-exp[-(x + x x - 2x D ) 2 /3kt]}dx 1 . (A29) 

and 

X(R,t) = v{R,t)R-W + i: A (J^) , (A30) 

where the first term is the disk evolution with a zero boundary condition and the second term is the fixed E boundary 
effect. As t4 oo, E -> E A (R/R D )- 3 / 2 . 

Figure 15 shows the evolution of the constant a = 0.01 disks with two different boundary conditions: E(10 AU)=50 
g cm~ 2 (black curves) and E(10 AU)=0 g cm -2 (red curves). The initial conditions are set as 

E(i?< 50AU,t = 0) =50gcm- 2 , (A31) 

E(i? > 50AU,t = 0) = OgcuT 2 . (A32) 

As shown in figure 15, the disk with zero boundary condition behaves quite similar to the similarity solution. However, 
for the £(10 AU)=50 g cm -2 disk, the influence of the boundary becomes significant with a part E oc R~ 15 and a 
outer part decreasing exponentially. 
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The mass accretion rate at Rd decreases to zero eventually (E oc R n 1 / 2 ). This can be shown by assuming 
E = gR m and v — kR n . Inserting these into 

M = 6irR 1/2 -^-(isZR 1/2 ) , (A33) 
oR 



we find 



M = 6nkg(m + n+ ^)R m+n . (A34) 



In the asymptotic case as discussed above with E oc R n 1 / 2 , M — 0, which means at t — > oo, M — > at the inner 
boundary R^, so no mass flow in the disk. This can be simply understood because the disk has zero torque with 
EociT"- 1 / 2 if v = kR n . 

However, notice that this E oc i?~ 3 / 2 behavior is only observed for layered disks with small or negligible dead zone 
viscosity, so that at very late stage Rd is very small and almost fixed. 

A. 3 Residual dead zone viscosity 
If the residual dead zone viscosity a r d >10 -4 , the dead zone can transport the mass at a rate comparable to the 
accretion rate of the active layer, considering Mi oc EjCtj and E,j >> E a . If we assume the disk has a constant mass 
accretion rate after a long period of time, then M—-2ir'RT la vii a -2'KTiEdVRd is a constant and we can integrate the 
angular momentum equation to derive 

(E a v Ra + E d v Rd )R 3 n = (i/„£ a + ^E d )i? 3 ^ + C . (A35) 

If we assume at the stellar surface i?* dVl/dR is 0, we find C—-MVL R 2 /27r. Thus 



M 

v a E a + VdX*d = — 
3tt 



R n 1/2 



(A36) 



Furthermore, if the disk is irradiation-dominated so that the temperature is ocR 1 / 2 , we have ^ a =k a R and ^^k^R 
where k a and k^ are constants and proportional to a a and ad- So 



^ M 



inkdR 



R n 1/2 



^E a . (A37) 
kd 



If the dead zone transports mass at a rate higher than the active layer, the last term can be neglected for order of 
magnitude argument. And compared with a constant-a disk having the same mass accretion rate, the layered disk has 
a surface density higher than the constant-a disk by a factor of a /ad- If we assume the outer pure MRI active disk 
and the inner layered disk has the same mass accretion rate, and the outer disk only has v a while the inner dead zone 
is dominated by the disk's surface density will increase by a factor of a a /ad within the dead zone outer radius R_d, 
which is shown in Figure 9. Also this density change at R_d is more gradual (equation 22 is a smooth function of R) 
than the zero dead zone viscosity case. 

B. GI DISK FRAGMENTATION RADIUS 

Gammie (2001) has pointed out that when the disk cooling timescale t coo i < 3f2 -1 the disk will fragment. By 
assuming local dissipation, Gammie (2001) has shown that 

4 

tco ° l = 9 7 (7-l)afi ' (B1) 

Thus, the disk will fragment if a > 0.06 (Rice et al. 2005). However, the above fragmentation condition has only 
been tested for disks without any irradiation. With irradiation dominated disk, t coo i is hard to be defined and a >1 
condition are used by Rafikov (2009) instead (Here we use a >0.06 to be consistent with the non-irradiated case.). 
For a constant accretion rate we have 

a n* = ^- (B2) 



Combined with Q=1.5, we derive 



a = ~2c!~ ' (B3) 



Thus the disk will fragment if MG/2c 3 >0.06. 
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For a viscous heating dominated Q=l disk, at a given mass accretion rate (M), the relationship between T c and R 
is given in Equation (23) in Paper I. If we reorganize the equation and assume (3—0 (k=CT q P^) with Q=3/2, we got 

2/(7-2a) 



ji _g2/(7-2a)2-12/(7-2a)^-9/(7-2a) 



M 



c 2/(7-2 a ) 



l/(7-2a) 



G ,l/(7-2a) M 3/(7-2a) 7r -4/(7-2a) ^ 



(B4) 



Inserting this into MG/2cf>0.06, 



i?>0.5 x .12( 14 - 4 «)/ 27 7r- 4 / 9 G(- 11+4 «)/ 27 C 2 / 9 a- 2 / 9 

k s(8-2a)/9 



M 1/3 M ( - 8+4a)/27 . 



(B5) 

Thus the fragmen- 



With our dust opacity a=0.738 (Zhu et al. 2009 b), the critical fragmentation radius RocM ° 19 . 

tation radius is insensitive to M in the viscous heating dominated case and ^100 AU (Figure 5). 

For an irradiation dominated case (low mass accretion rate), c s in equation B3 is determined by the irradiation 
(equation 4), thus 

-4/3 / k \ 2 / fL \ 1/2 

(B6) 



R > 0.06 ( MG) 



J \ fj,m,h J V 47rcr 



In this case the critical fragmentation R has a sharper dependence on M 4 / 3 (Figure 5). The GI fragmentation region 
is outlined in Figure 5, which agrees with Figure 10 of Cossins et al. (2009)b. 
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Fig. 1. — The disk mass accretion rate (upper panel) and the mass of the central star and disk (lower panel) with time for our fiducial 
model (Q c = l(T 14 rad s" 1 ). The upper curve in the lower panel represents the central star mass while the lower one represents the disk 
mass. 
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Fig. 2. — The disk surface density, temperature structure, mass accretion rate and a and Q-parameter at the end of the infall phase (0.2 
Myr) for our fiducial model ( solid curves, Sa = lOOg cm -2 and qm=0.01) and the model with a less massive active layer but a stronger 
MRI viscosity (dotted curves, = lOg cm -2 and om=0.1). The triangle labels the centrifugal radius at this time. In the upper left panel, 
the curves with higher surface densities within 20 AU represent the dead zone surface density, while the ones with E <100 g cm -2 and 
extend much further represent the active layer surface density. The upper curves in the upper right panel show the midplane temperature 
while the lower curves show the active layer temperature. The black curve in the lower left panel shows the disk mass accretion rate, while 
the red curve shows the portion of the mass accretion rate lead by the infall (the third term in equation 3). The upper curve in the lower 
right panel shows the Q parameter while the lower curve shows the midplane a (caused by GI). 
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Fig. 4. — The disk mass accretion rates and the mass of the central star and disk with time for different core rotations. The upper two 
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Fig. 5. — Unstable regions in the R — M plane for an accretion disk around IMq central star. The shaded region in the upper left shows 
the region subject to classical thermal instability. The shaded region in the middle shows where the central temperature of steady GI 
models exceeds an assumed MRI trigger temperature of 1400 K, thus subject to MRI-GI instability. The lines labeled Rm and Rq are the 
maximum MRI and minimum GI steady accretion radius (the boundaries of the shaded region, detail in Zhu et al. 2009b). The shaded 
region on the right shows where the GI disk will fragment (Sec Appendix for details). The fragmentation limit agrees with Cossins et al. 
(2009)b. 
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Fig. 7. — The disk surface density distributions at 1 million years for layered models (solid curves) and constant-o models (o=0.01, 
dotted curves) for f! = 2 X 10 — 14 rad s — 1 (red curves), fl = 10~ 14 rad s — 1 (black curves), and fl = 3 X 10 _15 rad s -1 (blue curves). The 
two black lines indicate S oc R~ 1 and S oc R~ 1,5 respectively. The right panels show the disks' mass beyond R (M(R)= J^°27rR£dR). 
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Fig. 8. — Disk evolution of the model without the active layer, in which case only GI has been considered to transport mass from the 
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rates and the upper right panel shows the masses of the central star and disk with time. The lower left panel shows the disk surface density 
distribution at 0.3 million years (solid curve) and 1 million years (dashed curve). The lower right panel shows the disk a and Q parameter 
at 0.3 and 1 Myr. 
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Fig. 9. — Disk evolution of the model similar to our fiducial layered model but with a residual dead zone viscosity (a r d = The 
upper left panel shows the disk mass accretion rates and the upper right panel shows the masses of the central star and disk with time. 
The lower left panel shows the disk surface density at 0.3 million years (solid curve) and 1 million years (dashed curve). The lower right 
panel shows the disk a and Q parameter at 0.3 and 1 million years. 
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R [AU] R [AU] 

Fig. 10. — For constant-^ disk (long dashed curve a= lCT 3 ; solid curve q=0.01; dotted curve a=0.1) with f! c =2xl0 14 rotating core, 
upper panels show the disk mass accretion rate (upper left panel) and central star mass (upper right panel) with time, while lower panels 
show the disk surface densities at 0.2 million years (lower left panel) and 1 million years (lower right panel) . The dash-dotted line in the 
lower left panel shows the surface density corresponding to Q=l disk. 
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Fig. 11. — In the layered disk model, the dead zone outer radius Rc with time for fl c = 2 X 10 14 rad s 1 (upper two curves) 
and Sl c = 3 X 10~ 15 rad s _1 (lower two curves). Two different configurations of the layered disk have been considered (Ea&m=^ and 
E a c*m=0.1). 
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Fig. 12. — Our simulated disk surface density compared with different MMSN models. The black solid curve is the surface density from 
one of our layered models (Q = 2 X 10 _14 rad s _1 at 0.6 Myr) compared with the surface density of the primordial disk constrained by 
Desch (2007). The red dashed curve is the surface density from our S7 = 10 _14 rad s _1 model at 2 Myr compared with the Minimum mass 
solar nebulae from Weidenschilling (1977) (dotted red curve). 
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Fig. 13. — The disk surface density distribution of our fiducial model (solid curve), fiducial model with a residual dead zone viscosity 
(a r( i = 10~ 3 , dashed curve), a constant-a (a=0.01) model (dotted curve) and a pure-GI model(long dashed curve) at 1 Myr compared 
with the disk surface densities constrained by millimeter observations of Ophiuchus (colored curves, Andrews et al. 2009). 
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Fig. 14. — The disk surface brightness temperature for the models in Figure 13: the fiducial model (solid curve), the fiducial model with 
a dead zone residual viscosity (a r[ i = 10 -3 ) (dashed curve), and a constant-o model with o=0.01 (dotted curve) at 1 millimeter (black 
curves) and 10 millimeters (red curves). EVLA will be able distinguish these different disk structures at 10 millimeters. 
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Fig. 15. — The viscous disk surface density with radii at t=50, 500, 5000 kyr (curves from up to bottom) with the fixed S inner boundary 
condition (Black curves for the inner boundary S=50 g cm -2 while the red curves for S=0 g cm -2 ). The dotted curve extending to 50 AU 
shows the disk initial condition, while the dashed line shows the asymptotic limit p ocR — 1 ' 5 . The open dots at the top shows the scaling 
radius (see the definition in Lynden-Bell & Pringle 1974) Rt+10 AU at t=50, 500, 5000 kyr. 
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TABLE 1 

1D2Z MODELS AT 0.3 AND 1 MYR 



cloud fl 




OL M 


Rd 0.3/1 Myr 


M ln a 


M b 


M t c 






g cm -2 




AU 


M 


M s 


M Q 


2x10" 


■14 


100 


0.01 


40/13 


0.29/0.19 


0.18/0.14 


0.47/0.33 


2xl0~ 


■14 


10 


0.01 


127/66 


0.42/0.33 


0.05/0.09 


0.48/0.42 


2x10" 


■14 





0.01 


-/- 


0.47/0.41 


-/- 


0.47/0.41 


1x10" 


■14 


100 


0.01 


24/7.6 


0.17/0.1 


0.07/0.06 


0.24/0.16 


1x10" 


■14 


10 


0.01 


47/34 


0.32/0.23 


0.008/0.022 


0.32/0.26 


3xl0~ 


■15 


100 


0.01 


4.6/2.4 


0.05/0.04 


0.006/0.007 


0.057/0.049 


3x10" 


■15 


10 


0.01 


5.5/4.6 


0.06/0.06 


0.0004/0.001 


0.063/0.063 


3x10- 


15 





0.01 


-/- 


0.064/0.064 


-/- 


0.064/0.064 



a The disk mass within 
b The disk mass beyond Rd 
c Thc total disk mass 



